{
 "cells": [
  {
   "cell_type": "markdown",
   "id": "816baa7b",
   "metadata": {
    "lines_to_next_cell": 0
   },
   "source": [
    "#  9\\.  回归方法综合应用练习  # \n",
    "\n",
    "##  9.1.  介绍  # \n",
    "\n",
    "本次挑战将会结合前面学习到的相关回归分析方法，完成一个多元回归分析任务。同时，我们将结合 NumPy，scikit-learn，SciPy，statsmodels 等库，复习不同方法的实现及应用。 \n",
    "\n",
    "##  9.2.  知识点  # \n",
    "\n",
    "  * 一元线性回归 \n",
    "\n",
    "  * 多元线性回归 \n",
    "\n",
    "  * 假设检验 \n",
    "\n",
    "挑战使用《An Introduction to Statistical Learning》书中提供的 ` Advertising  ` 示例数据集用于练习。该书是统计学习经典著作，有兴趣的话可以从作者网站 [ 免费下载阅读 ](http://www-bcf.usc.edu/~gareth/ISL/) 。 \n",
    "\n",
    "首先，我们加载并预览该数据集。 "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "ea5490e0",
   "metadata": {
    "lines_to_next_cell": 0
   },
   "outputs": [],
   "source": [
    "# 下载数据集\n",
    "!wget -nc https://cdn.aibydoing.com/aibydoing/files/advertising.csv"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "44070cbf",
   "metadata": {},
   "outputs": [],
   "source": [
    "import pandas as pd\n",
    "\n",
    "data = pd.read_csv(\"advertising.csv\", index_col=0,)\n",
    "data.head()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "0461064b",
   "metadata": {
    "lines_to_next_cell": 0
   },
   "source": [
    "|  tv  |  radio  |  newspaper  |  sales   \n",
    "---|---|---|---|---  \n",
    "1  |  230.1  |  37.8  |  69.2  |  22.1   \n",
    "2  |  44.5  |  39.3  |  45.1  |  10.4   \n",
    "3  |  17.2  |  45.9  |  69.3  |  9.3   \n",
    "4  |  151.5  |  41.3  |  58.5  |  18.5   \n",
    "5  |  180.8  |  10.8  |  58.4  |  12.9   \n",
    "  \n",
    "数据集包含 4 列，共 200 行。每个样本代表某超市销售相应单位件商品所需要支出的广告费用。以第一行为例，表示该超市平均销售 22.1 件商品，需要支出的电视广告费用，广播广告费用以及报刊广告费用为：230.1 美元，37.8 美元和 69.2 美元。 \n",
    "\n",
    "所以，本次挑战将前 3 列视作特征，最后一列视作目标值。 \n",
    "\n",
    "Exercise 9.1 \n",
    "\n",
    "挑战：使用 SciPy 提供的普通最小二乘法分别计算 3 个特征与目标之间的一元线性回归模型拟合参数。 \n",
    "\n",
    "规定：需使用 ` scipy.optimize.leastsq  ` 函数完成计算，并直接输出函数结果无需处理。 "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "8ef27c84",
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "from scipy.optimize import leastsq\n",
    "\n",
    "## 代码开始 ### (≈ 10 行代码)\n",
    "params_tv = None\n",
    "params_radio = None\n",
    "params_newspaper = None\n",
    "## 代码结束 ###\n",
    "\n",
    "params_tv[0], params_radio[0], params_newspaper[0]"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "27f7240a",
   "metadata": {
    "lines_to_next_cell": 0
   },
   "source": [
    "参考答案  Exercise 9.1 "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "221cc0c2",
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "from scipy.optimize import leastsq\n",
    "\n",
    "### 代码开始 ### (≈ 10 行代码)\n",
    "p_init = np.random.randn(2)\n",
    "\n",
    "def func(p, x):\n",
    "    w0, w1 = p\n",
    "    f = w0 + w1*x\n",
    "    return f\n",
    "\n",
    "def err_func(p, x, y):\n",
    "    ret = func(p, x) - y\n",
    "    return ret\n",
    "\n",
    "params_tv = leastsq(err_func, p_init, args=(data.tv, data.sales))\n",
    "params_radio = leastsq(err_func, p_init, args=(data.radio, data.sales))\n",
    "params_newspaper = leastsq(err_func, p_init, args=(data.newspaper, data.sales))\n",
    "### 代码结束 ###\n",
    "\n",
    "params_tv[0], params_radio[0], params_newspaper[0]"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "90322750",
   "metadata": {
    "lines_to_next_cell": 0
   },
   "source": [
    "期望输出 "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "0b5d19ce",
   "metadata": {},
   "outputs": [],
   "source": [
    "((array([7.03259354, 0.04753664]),\n",
    " (array([9.31163807, 0.20249578]),\n",
    " (array([12.35140711, 0.0546931]))"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "dfb7cbf6",
   "metadata": {
    "lines_to_next_cell": 0
   },
   "source": [
    "接下来，我们根据最小二乘法求得的结果，将拟合直线绘制到原分布散点图中。 \n",
    "\n",
    "Exercise 9.2 \n",
    "\n",
    "挑战：以横向子图的方式绘制 3 个特征分别与目标之间的散点图，并添加线性拟合直线。 \n",
    "\n",
    "规定：线性拟合直线开始于散点图中最小横坐标值，结束于最大横坐标值，并以红色显示。 "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "02d153d1",
   "metadata": {},
   "outputs": [],
   "source": [
    "from matplotlib import pyplot as plt\n",
    "\n",
    "%matplotlib inline\n",
    "\n",
    "## 代码开始 ### (≈ 10 行代码)\n",
    "\n",
    "## 代码结束 ###"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "e11218b3",
   "metadata": {
    "lines_to_next_cell": 0
   },
   "source": [
    "参考答案  Exercise 9.2 "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "0c46fc56",
   "metadata": {},
   "outputs": [],
   "source": [
    "from matplotlib import pyplot as plt\n",
    "%matplotlib inline\n",
    "\n",
    "### 代码开始 ### (≈ 10 行代码)\n",
    "fig, axes = plt.subplots(1, 3, figsize=(16, 4))\n",
    "\n",
    "data.plot(kind='scatter', x='tv', y='sales', ax=axes[0])\n",
    "data.plot(kind='scatter', x='radio', y='sales', ax=axes[1])\n",
    "data.plot(kind='scatter', x='newspaper', y='sales', ax=axes[2])\n",
    "\n",
    "x_tv = np.array([data.tv.min(), data.tv.max()])\n",
    "axes[0].plot(x_tv, params_tv[0][1]*x_tv + params_tv[0][0], 'r')\n",
    "\n",
    "x_radio = np.array([data.radio.min(), data.radio.max()])\n",
    "axes[1].plot(x_radio, params_radio[0][1]*x_radio + params_radio[0][0], 'r')\n",
    "\n",
    "x_newspaper = np.array([data.newspaper.min(), data.newspaper.max()])\n",
    "axes[2].plot(x_newspaper, params_newspaper[0][1] *\n",
    "             x_newspaper + params_newspaper[0][0], 'r')\n",
    "### 代码结束 ###"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "beb9fe53",
   "metadata": {
    "lines_to_next_cell": 0
   },
   "source": [
    "期望输出 \n",
    "\n",
    "![image](https://cdn.aibydoing.com/aibydoing/images/document-uid214893labid7506timestamp1546077405080.png)\n",
    "\n",
    "接下来，我们尝试建立包含全部特征的多元线性回归模型。 \n",
    "\n",
    "$$ y = \\omega_0 + \\omega_1 * tv + \\omega_2 * radio + \\omega_3 * newspaper $$ \n",
    "\n",
    "Exercise 9.3 \n",
    "\n",
    "挑战：使用 scikit-learn 提供的线性回归方法建立由 3 个特征与目标组成的多元线性回归模型。 \n",
    "\n",
    "规定：仅能使用 scikit-learn 提供的 ` sklearn.linear_model.LinearRegression  ` 类。 "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "6599032d",
   "metadata": {},
   "outputs": [],
   "source": [
    "from sklearn.linear_model import LinearRegression\n",
    "\n",
    "## 代码开始 ### (≈ 4 行代码)\n",
    "model = None\n",
    "## 代码结束 ###\n",
    "\n",
    "model.coef_, model.intercept_  # 返回模型自变量系数和截距项"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "ac4bf094",
   "metadata": {
    "lines_to_next_cell": 0
   },
   "source": [
    "参考答案  Exercise 9.3 "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "c4c14fd4",
   "metadata": {},
   "outputs": [],
   "source": [
    "from sklearn.linear_model import LinearRegression\n",
    "\n",
    "### 代码开始 ### (≈ 4 行代码)\n",
    "X = data[['tv', 'radio', 'newspaper']]\n",
    "y = data.sales\n",
    "\n",
    "model = LinearRegression()\n",
    "model.fit(X, y)\n",
    "### 代码结束 ###\n",
    "\n",
    "model.coef_, model.intercept_  # 返回模型自变量系数和截距项"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "a5a5d0de",
   "metadata": {
    "lines_to_next_cell": 0
   },
   "source": [
    "期望输出 "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "1b03cddd",
   "metadata": {},
   "outputs": [],
   "source": [
    "(array([ 0.04576465,  0.18853002, -0.00103749]), 2.9388893694594067)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "e05948bb",
   "metadata": {
    "lines_to_next_cell": 0
   },
   "source": [
    "接下来，我们希望对多元线性回归模型进行检验。使用 statsmodels 库提供的相关方法来完成拟合优度检验和变量显著性检验。 \n",
    "\n",
    "Exercise 9.4 \n",
    "\n",
    "挑战：使用 statsmodels 库提供的相关方法来完成上面多元回归模型的拟合优度检验和变量显著性检验。 \n",
    "\n",
    "提示：你可以使用 ` statsmodels.api.sm.OLS  ` 或 ` statsmodels.formula.api.smf  ` ，后一种实验未涉及，需自行了解学习。 "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "eb61fd80",
   "metadata": {},
   "outputs": [],
   "source": [
    "import statsmodels.formula.api as smf\n",
    "\n",
    "## 代码开始 ### (≈ 3 行代码)\n",
    "results = None\n",
    "## 代码结束 ###\n",
    "\n",
    "results.summary2()  # 输出模型摘要"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "cd7a6b10",
   "metadata": {
    "lines_to_next_cell": 0
   },
   "source": [
    "参考答案  Exercise 9.4 "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "2b84ea96",
   "metadata": {},
   "outputs": [],
   "source": [
    "import statsmodels.formula.api as smf\n",
    "\n",
    "### 代码开始 ### (≈ 3 行代码)\n",
    "results = smf.ols(formula='sales ~ tv + radio + newspaper', data=data).fit()\n",
    "### 代码结束 ###\n",
    "\n",
    "results.summary2() # 输出模型摘要"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "f13ad5a7",
   "metadata": {},
   "source": [
    "期望输出 \n",
    "\n",
    "Model:  |  OLS  |  Adj. R-squared:  |  0.896   \n",
    "---|---|---|---  \n",
    "Dependent Variable:  |  sales  |  AIC:  |  780.3622   \n",
    "Date:  |  |  BIC:  |  793.5555   \n",
    "No. Observations:  |  200  |  Log-Likelihood:  |  -386.18   \n",
    "Df Model:  |  3  |  F-statistic:  |  570.3   \n",
    "Df Residuals:  |  196  |  Prob (F-statistic):  |  1.58e-96   \n",
    "R-squared:  |  0.897  |  Scale:  |  2.8409   \n",
    "|  Coef.  |  Std.Err.  |  t  |  P>|t|  |  [0.025  |  0.975]   \n",
    "---|---|---|---|---|---|---  \n",
    "Intercept  |  2.9389  |  0.3119  |  9.4223  |  0.0000  |  2.3238  |  3.5540   \n",
    "tv  |  0.0458  |  0.0014  |  32.8086  |  0.0000  |  0.0430  |  0.0485   \n",
    "radio  |  0.1885  |  0.0086  |  21.8935  |  0.0000  |  0.1715  |  0.2055   \n",
    "newspaper  |  -0.0010  |  0.0059  |  -0.1767  |  0.8599  |  -0.0126  |  0.0105   \n",
    "Omnibus:  |  60.414  |  Durbin-Watson:  |  2.084   \n",
    "---|---|---|---  \n",
    "Prob(Omnibus):  |  0.000  |  Jarque-Bera (JB):  |  151.241   \n",
    "Skew:  |  -1.327  |  Prob(JB):  |  0.000   \n",
    "Kurtosis:  |  6.332  |  Condition No.:  |  454   \n",
    "  \n",
    "我们可以看到，这里得到的回归拟合系数和上文 scikit-learn 计算结果一致。于此同时，tv 和 radio 的 P 值接近于 0 [精度]，而 newspaper 的 P 值相对较大。我们可以认为 tv 和 radio 通过了变量显著性检验，而 newspaper 则未通过。实际上，你可以尝试去掉 newspaper 特征之后重新计算多元线性回归结果。 \n",
    "\n",
    "Model:  |  OLS  |  Adj. R-squared:  |  0.896   \n",
    "---|---|---|---  \n",
    "Dependent Variable:  |  sales  |  AIC:  |  778.3941   \n",
    "Date:  |  |  BIC:  |  788.2891   \n",
    "No. Observations:  |  200  |  Log-Likelihood:  |  -386.20   \n",
    "Df Model:  |  2  |  F-statistic:  |  859.6   \n",
    "Df Residuals:  |  197  |  Prob (F-statistic):  |  4.83e-98   \n",
    "R-squared:  |  0.897  |  Scale:  |  2.8270   \n",
    "|  Coef.  |  Std.Err.  |  t  |  P>|t|  |  [0.025  |  0.975]   \n",
    "---|---|---|---|---|---|---  \n",
    "Intercept  |  2.9211  |  0.2945  |  9.9192  |  0.0000  |  2.3403  |  3.5019   \n",
    "tv  |  0.0458  |  0.0014  |  32.9087  |  0.0000  |  0.0430  |  0.0485   \n",
    "radio  |  0.1880  |  0.0080  |  23.3824  |  0.0000  |  0.1721  |  0.2038   \n",
    "Omnibus:  |  60.022  |  Durbin-Watson:  |  2.081   \n",
    "---|---|---|---  \n",
    "Prob(Omnibus):  |  0.000  |  Jarque-Bera (JB):  |  148.679   \n",
    "Skew:  |  -1.323  |  Prob(JB):  |  0.000   \n",
    "Kurtosis:  |  6.292  |  Condition No.:  |  425   \n",
    "  \n",
    "上面是去掉 newspaper 特征之后重新计算的结果，你可以发现包含和未包含 newspaper 的多元线性回归模型  $R^2$  值均为 0.896。这也就印证了该特征对于反映目标数值变化并无太大帮助。你可以尝试去除其余两个特征之一后，再查看  $R^2$  拟合优度值的变化。 "
   ]
  }
 ],
 "metadata": {
  "jupytext": {
   "cell_metadata_filter": "-all",
   "main_language": "python",
   "notebook_metadata_filter": "-all"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}
